Iranian Journal of Numerical Analysis and Optimization onze 


Vol. 9, No. 1, (2019), pp 1-16 5 
DOI:10.22067 /ijnao.v9i1.53910 


Transformation to a fixed domain in 
LP modelling for a class of optimal 
shape design problems 
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Abstract 


A class of optimal shape design problems is studied in this paper. The 
boundary shape of a domain is determined such that the solution of the 
underlying partial differential equation matches, as well as possible, a given 
desired state. In the original optimal shape design problem, the variable 
domain is parameterized by a class of functions in such a way that the optimal 
design problem is changed to an optimal control problem on a fixed domain. 
Then, the resulting distributed control problem is embedded in a measure 
theoretical form, in fact, an infinite-dimensional linear programming problem. 
The optimal measure representing the optimal shape is approximated by a 
solution of a finite-dimensional linear programming problem. The method is 
evaluated via a numerical example. 


Keywords: Approximation; Optimal shape design; Linear programming; 
Measure theory. 


1 Introduction 


An optimal shape design (OSD) problem is concerned with the optimization 
of a performance index depending on the shape of some region. Many-faceted 
problems naturally arise in enginecring applications with the goal of designing 
a specific structure in an optimal sense. Typical applications are design of a 
nozzle [12], airfoil boundary [4,10], and spacecraft shape [6] with respect to 
specific optimality conditions. 
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There are several methods for a numerical solution of an OSD problem, 
for example, gradient based methods [2], the Newton method [9,15], sequen- 
tial quadratic method [19], and fictitious domain method [8]. In the above 
mentioned methods, computation of solution is a time-consuming task, be- 
cause they need to solve many boundary value problems. Moreover, these 
methods are iterative and they need to have an initial guess for solution. 

The main focus of this article is to find an appropriate formulation of 
the optimal shape design problem that is attractive for consistent numerical 
computation. The advantages of the proposed method lie in the fact that the 
method is not iterative, that it is self-starting, and that it does not need to 
solve the corresponding boundary value problem. Moreover, in the model- 
ing of the OSD problem, differentiability of the cost function being a limiting 
condition in the above mentioned methods, is reduced to a measurability con- 
dition. In the methods of [9], the control and state functions are discretized 
and the linearity of variational forms is an inherent condition, while it is not 
a restrictive condition in the proposed method, where a general variational 
form is studied. 

In the present work, the variable domain is first parameterized by a class of 
functions in such a way that the OSD problem changes to an optimal control 
problem accompanied with partial differential equations on a fixed domain. 
The converted problem can be solved by extending the measure theoretical 
method proposed by Rubio [16]. The methods of this type were frequently 
and successfully applied to different OSD problems. In [3] the problem of 
finding the upper wall of a nozzle was discussed where the goal is to match 
the velocity of flow to a constant prescribed value in a given subregion near 
exit. It also includes details about the measure theoretical approach to solve 
OSD problems. The same problem for pressure of the flow was given in [12] 
with more numerical simulations to study the sensitivity analysis and stability 
of the method. In [4] the problem of finding the optimum shape for a low 
speed airfoil has been solved by the method. A problem concerning with the 
control of thermoelastic deformation was also presented in [13], where the 
extension of the method of [3] has been used to solve it. 

The main difference of the present work with the above mentioned re- 
searches is that the optimization occurs in a completely unknown two- 
dimensional area and that there is no restriction to any special form of the 
problem. Such a generalization of the domain is also presented in [7], which 
is based on the domain discretization, while in the present paper, the trans- 
formation of domain is used. 


2 The optimal shape design problem 


The set of admissible shapes in the definition of an OSD problem is the set of 
all possible shapes in which we seek for the best. In our problem definition, let 
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consider a class of OSD problems where the admissible shapes are described 
by two-dimensional domains with a free boundary: 


O = {& = (%1,%2)7 CR? | 2 EI, #1 € (0, u(%2))}, 


where J = (0,1). The set of admissible shapes is usually parameterized or 
discretized in order to reduce the search process from the class of admissi- 
ble shapes to a set of admissible functions. So, we parameterize the above 
domains by a function u € Ugg, where 


2 GN laxaaly, 


Uaa = {u € C*(I) | u(0) = u(1) = 1, u(%2) € [Ar, Ba], is 


where £3; > 0, 82, a1, and a are given. 


By measuring the desirability of admissible shapes, let J be a cost function 
depending on a function y and its gradient as follows: 


i= [ favanae, (1) 


where f is a nonlinear real-valued measurable function and y is defined by 
the solution of a second order elliptic boundary value problem on 9. The 
boundaries of the feasible domains are divided into the moving boundary 
part 

I™ = {(41,%2)7 € R? | &, = u(%e) for all % € J}, 


and the fixed boundary part I” with 02 =P” UT and IM ATF = 0. 


Now, the OSD problem is defined by 


in J 
sas? 
with the equality constraint 
| a(¥, n)dz = | i(q)dT for all 9 € V, 
Q an 


which is a variational form of the underlying PDE, V Cc H 1Q) is the space 
of corresponding test functions, and H!(Q) is a Sobolev space of square 
integrable functions with square integrable first derivatives on Q. Functions 
@ and / are assumed to be nonlinear in general. We use the notations of [9] 
in our article. 
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Xj 


Figure 1: Transformation of the domain. 


3 Transformation to a fixed domain 


To cope with the variation of the domain, we transform the optimal shape 
design problem into an optimal control problem on a fixed domain. This 
routine approach has been used in literature before, for example [9, 13]. 

As sketched in Figure 1, the transformation defined by 


T3203) ; ; 
(#1, %2)¢ + (1,22) = (=) (2) 


maps the moving domain {2 to the fixed reference area 2. = I x I. 
So, we consider the problem over a fixed domain 


in J(u, 
join Jy), 
with the equality constraint 
du 
a(y,7, U, — )dx = ly for all ne V, (3) 
Q dx 


where 1, is a function depending only on the test function 7. Here, a includes 
all terms depending on y, 7, u, and . 

By this transformation, the unknown geometrical factors are transformed 
from boundary to equations, that is, the equations must be solved on a fixed 
known domain while they contain unknown parameters. It will be shown in 
the next sections that this type of problems may be solved by the described 
method. 

The transformation will also changes the integrand of functional (1) as a 
function f dependent on y, Vy, u and ee: 
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du 
J= if f(u,4 VY, Ce a (4) 


Now, we choose the derivative of u(x2) as a control variable w in order to 
obtain the following control system: 


du 


=— WwW 
dx , 


u(0) = u(1) = 1. (6) 


ae. 2 ET, (5) 


The differential form of this control system should be changed as a variational 
form consistent to (3). To do this, let define an admissible pair for (5)—(6). 


Definition 1. The pair p = (w,u) is said admissible if the control function 
w:I — [a1,a2] is measurable and u € Usa is the response of the system 
(5)-(6) to w. 


As unknown factors in (3)-(4) are depending to u, then for a given admis- 
sible pair p = (w,u), there exists a solution for (3) with a cost measured by 
(4). So, the problem will be finding optimal admissible pair in the set of all 
admissible pairs which is assumed nonempty and is denoted by P. However 
to find the best admissible pair we have to determine characteristic properties 
of an admissible pair. 

Let TY = 2 x [(81, Be] x [a1,a2] x A, where A C R? is the set that 
(y(#1, %2), Vy(%1,22)) gets its values in it. Moreover, assume that B is an 
open ball in R? containing [(1, 82] x I, and let C1(B) be the space of all 
real-valued continuous differentiable functions with continuous first partial 
derivatives on B. For y € C1(B) and p = (w,u) € P, define yp” € C(Y) as 
follows: 


yp" (£1, £2, u(X2), w(X2),y, Vy) = pi(u(x2), 2)w(e2) + ya(u(x2), 22). (7) 
In (7), Y1, %2 denote the first partial derivatives of y with respect to its 
variables. Admissibility of p = (w, u) implies: 


1 
| yp (#1, £2, u(2), w(r2),y, Vy)dx2 
0 


7 | {pr (uaa), 22)w(a2) + o2(u(2), #2) }derg 
Soul GOO) Ae. 


for all y € C'(B), where y(u(1),1) and y(u(0),0) are known. Integrating 
the above equation with respect to x2 on J, the integral form of (5)—(6) can 
be written as 


7. py’ (x1, £2, u(£2), w(r2),y, Vy)de = Ay for ally €C*(B). (8) 
) 
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Therefore, the optimal control formulation of the OSD problem may be in- 
terpreted as the minimization of (4) over variational constraints (3) and (8). 


4 Measure theoretical approach 


By Definition 1, each shape corresponds to exactly one admissible pair p = 
(w,u) € P. This is an injection between U,g and P. So, the minimization 
of (1) over the class of shapes is equivalent to the minimization of (4) over 
P. Now, we transfer the problem of minimization of (4) over P into another 
nonclassical problem. Let us define the following mapping: 


hpi Fo | re, woiieiavem, Pec, 
Q 


where C(Y) is the space of all continuous real-valued functions defined on 
Y. As an injection, the transformation p ++ A, provides us to describe 
the set of all admissible pairs P as a subset of the set of all linear continuous 
mappings on C(Y). Moreover, by the Riesz representation theorem (see [18]), 
corresponding to Ap, there is a positive Borel measure y on Y such that 


A,(F) = w(F), for all F € O(Y). 


Now, the problem of minimization of the functional (4) over P is enlarged 
to the minimization of 


LT: wwf) (9) 

over the set of all measures yz satisfying 
ple) tay for all n € V, (10) 
U(Gy) = Ay, for all y € C'(B), (11) 


where F,,(#1, %2, u(x2), w(x), y, Vy) = aly, 7, u, a) and Gy = y™. In other 
words, (10) and (11) are a measure theoretical interpretation of (3) and (8) 
respectively. 


We denote the set of all positive Borel measures on YT satisfying (10)—(11) 
as Q. We also assume that M*(TY) is the set of all positive Borel measures 
on Y. If we consider the space M*(Y) with weak*-topology, it can be seen 
from [16], that Q is compact. In the sense of this topology, the functional 
I :Q— R defined by (9) is a linear continuous functional on the compact 
set Q, thus it attains its minimum on Q (see [16, Theorem III.1]), and then 
the measure-theoretical problem, which consists of finding the minimum of 
the functional (9) over a subset of M+(Y), possesses a minimizing solution, 
say p*, in Q. 
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The minimization problem of (9)—(11) is an infinite-dimensional LP prob- 
lem, and we are mainly interested in approximating it. It is possible to 
approximate the solution of the problem (9)—(11) by the solution of a finite- 
dimensional linear program of sufficiently large dimension. 

We first consider the minimization of (9) not only over the set Q but over 
a subset of it defined by requiring that only a finite number of constraints 
(10)—(11) are satisfied. 

Consider equalities (10)—(11), let the sets of functions {n;, i = 1,2,3,...} 
and {y;, j = 1,2,3,...} have dense linear combinations in V and C1(B), 
respectively. Now we can prove the following. 


Proposition 1. Let Q(M,, M2) be a subset of M*(Y) consisting of all mea- 
sures satisfying 


BF; ) les t= 1,2,..., My, 
(Ou Degg, GB ye: 


Tf 0(.M1, M2) = infgcm,,m,) H(f) and 0 = infg u(f), then 0(M1, M2) > 0 as 
My, M2 > ow. 


Proof. See the proof of [3, Proposition 2]. oO 


The first stage of approximation is completed successfully. As the second 
stage, it is possible to develop a finite-dimensional linear program whose 
solution can be used to construct the solution of the infinite-dimensional 
linear program of minimizing (9) subject to (10)—-(11). From [16, Theorem 
A.5], we can characterize a measure, say p*, in the set Q(M, M2) at which 
the function > p(f) attains its minimum. 


Proposition 2. The measure p* in the set Q(M,, M2) in which the function 
pu > Lf) attains its minimum has the form 


M,+Mz2 


w= S20 6(z2) (12) 


k=1 
with z, € T and of >0,k =1,2,...,M@i + Mo. 
Here 5(z) is the atomic measure concentrated at z, characterized by 


wensia{a $258 


for all S Cc Y; see [17]. This also implies that 6(z)(H) = H(z), where 
H € C(Y) and z € Y. Now, the measure theoretical optimization prob- 
lem is equivalent to a nonlinear optimization problem, in which unknowns 
are coefficients 97, and supports zj, fork = 1,2,..., Mi + Mz. It would be 
convenient if we could minimize the function w > p(f) only with respect 
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to the coefficients 97, fork = 1,2,...,My + Mg, in (12), which would be a 
linear programming problem. However, we do not know the support of the 
optimal measure. The answer lies in the approximation of this support by 
introducing a dense set in YT. 


Proposition 3. Let EF be a countable dense subset of T. Given € > 0, a 
measure \ € Mt(Y) can be found such that 


Iu" -— AA <e, 
I(e* — A) Fa, ) Se i= 1,2, Mi, 
[(u* — A)(Gy,)| Se J=41,4; ., Mp, 
and the measure X has the form 
Mi+M2 
A= DS) of (zn), 
k=1 


where the coefficients 0; are the same as in the optimal measure (12) and 
Zee E. 


Proof. See the proof of [16, Proposition II.3 ]. oO 


Thus, the infinite-dimensional linear programming (9) with restrictions 
defined by (10)-(11) can be approximated by another LP problem in which 
Ze, for k = 1,2,...,N, belongs to the known dense subset of TY. To 
construct such a subset, we grid entire of YT into N >> M, + Mp2 nodes 
z, = (rk, ok, u*,w*,y*, Vy"). This gird may be generated by dividing 0, 
(C1, G2], [a1, a2], and A separately and rearranging nodes from 1 to N. For 


example, one could divide [0,1] for 71 to ny points: 111, 712,...,21n, which 
repeat in arrangement of z,’s. 
Now, we have the following LP problem with 01, 02,..., an as nonnegative 


decision variables: 


N 
Minimize > of (Zk) 
k=1 


subject to 
N 
SS” op Fne tel = Ty, i=1,2,...,Mi, 
k=1 
N 
\ ewan: jg =1,2,...,Mo, 
k=1 


The solution of this LP problem determines nodes that make the optimal 
measure and finally the optimal shape. As the solution may be trivial (zero), 
it is necessary to prevent from trivial situation. To do this, we benefit from 
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Figure 2: A typical grid on 9. 


theoretical properties of measure. Let consider a simple example in which 
Q is covered with 48 nodes as given in Figure 2. If H is the characteristic 
function on the gray area in the picture, then the measure of H will be the 
area of band. On the other hand, this measure is equal to > 9%, where the 
summation is on the set of all points with «* equaling to 21;. So, we have a 
new constraint preventing the measure to be zero on this area: 


5 Ok = %12— 11 


{k:ak=211} 


This constraint forces the solution to be nonzero and preserving the properties 
of measure in the mentioned vertical band. We have to add similar constraints 
for all vertical and horizontal bands of this picture. 


4.1 Extraction of optimal control 


After solving the LP problem, we will find the optimal values for decision 
variables. The form of measure given in (12) is a combination of coefficients 
0, and support points z,. This tel us that, if a coefficient is zero, then the 
corresponding support does not have a role in the optimal measure. There- 
fore, nonzero variables will make the optimal measure in a fashion which we 
discuss here. In fact, we use the effect of this optimal measure to find the 
approximate solution for the control function. First, define yo = 0 and 


YL = S° Ok; 


keEA;, 
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m= Du 
|Ai| 


where 
A; = {k: on > 0, ak = 244}, iw ey eee One 


Now the approximation for optimal control function w(.) is a piecewise con- 
stant function that equals w; on each interval [y-1,y]. The trajectory is 
then simply found by solving the differential equation (5) with the initial 
condition u(0) = 1. The resulting solution is a piecewise linear function. 

We may also use numerical methods such as interpolation or fitting meth- 
ods to make a smooth shape. Choosing the convenient method depends on 
design purposes. A numerical treatment is also given in [12] which may be 
useful in this connection; see also [1] for numerical smoothing and approxi- 
mation methods. 


5 Numerical example 


In this section, a problem with applications in aerodynamic shape optimiza- 
tion is studied to validate the presented method from the computational 
standpoint. 


The velocity 0(%) at a point % € Q in a nonviscous incompressible poten- 
tial flow (such as for air or water at moderate speed) may be approximated 


by 7 
i@)=Vae), #4, 


where y satisfies the second order partial differential equation on Q: 
Ay = 0. (13) 


This type of potential flow occurs mainly in design of low speed airfoils, wings, 
and blades. The variational form of the above PDE is as follows (see [14]): 


| ViVidi = / AVy.ndl, for all 7 € V, (14) 
Qa rere) 
where 7 is the outward normal vector to T= 09. 


We also assume the following Dirichlet condition on the fixed boundary 
part 
ylize = 0. (15) 


So, the above variational form is changed to 


1 
i ViVidt = | Aydx., for all 7 EV, 
Q 0 
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where V = {n € H1(9): lar = 0}. 
Now, the specific transformation (2) can be used on the variable domain. 
This leads to the following transformed bilinear form 


1 atw? 
(-—+ )yim — ciwyem — ciwyine + uyene2| dx 
Joa 


1 
-| n(1, t2)yi (1, 22) V1 + wdao. 
Jo 


The above integral is now defined on the fixed domain 2, where coefficient 
functions depend on the parameter function u € U,q and its derivative w. 


It is easy to check that 7 = sinh(7%1) sin(7%2) is a solution of the bound- 
ary value problem (13), (15), or equivalently (14)—(15). 

To make a comparison between exact and numerical solution, we follow 
an inverse scheme to construct an OSD problem with known exact solution. 
To do this, we choose a cost function as 


iL Vg — Vola’, 
Tr 


which leads to a velocity matching problem. In other words, we want to 
find the moving boundary T™ in such a way that the velocity of the flow 
matches, as well as possible, a given velocity Vo. For applications of this type 
of matching OSD, we may address, for example, the pressure distribution 
matching in airfoil section design (see [5]). 

The above cost function finds the following format after transformation: 


1 
| |\Vy — Voll V1 + w2dae 
0 


or 


E (Sy = Valse. (16) 


In this example, for test functions of type 7, polynomials of x; and x2 with 
compact support on 2 of the following forms are chosen: 


ei(l—a2)*, #3(1-—21)%, Be De ses cct 


For test functions of type y, we use polynomials of u, functions depending 
only on a2, and functions in C!(B) with compact support on J; see [3, Sec.5] 
for more details about test functions. 

To illustrate the scheme, two specific moving boundaries are implemented: 

a: If the moving part of the boundary I™, is ¢, = 1, then the the above 
solution gives a velocity vector as 


Vy = (mcosh(zm) sin(7%2), 7 sinh(7) cos(7%2)) 
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Oo 
iho 
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Figure 3: Resulting control function, part a. 


on I™. Now, if we let Vo = (mcosh(z) sin(7%2), 7 sinh(m) cos(7#2)), then 
the OSD problem of minimizing (16) subject to (14)-(15) has the unique 
solution %; = 1 or equivalently u = 1. The solution of the corresponding 
linear programming problem is used to construct the control function w as 
shown in Figure 3. The corresponding optimal trajectory, u, is compared 
with the exact solution in Figure 4. It is clear that the numerical solution 
approximates the exact solution with a small error. 

In this case, we choose N = 30625, 6, = 0.4, 62 = 1.38, ay = —0.8, 
a2 = 0.9, My = 1, and Mz = 12. The LP problem has been solved by two 
phase revised simplex method, and the corresponding optimal objective was 
found as 0.01194. 

b: If the moving part of the boundary I™, is ¢; = 1+sin(maq), then the 
above solution gives a velocity distribution on I™ as 


Vy = (mcosh(m(1 + sin(ma#2))) sin(7%2), 7 sinh(m(1 + sin(a#2))) cos(7%2)). 
Now if we set 
Vo = (m cosh(r(1 + sin(7%2))) sin(w%2), 7 sinh(7(1 + sin(7%2))) cos(7¥2)), 


then the OSD problem of minimizing (16) over (14)—(15) has the unique solu- 
tion u = 14+sin(7%2). The solution of the corresponding linear programming 
problem is used to construct the control function w as shown in Figure 5. To 
enable a comparison, Figure 6 shows the resulting numerical shape and the 
exact solution of the problem that indicates the accuracy of results. 

In this part of the example, we chose N = 21875, 6, = 0.62, 62 = 1, 
a, = —1.96, ag = 1.96, M, = 1, and Mz = 12. The optimal objective of the 
corresponding LP problem was found as 2.1275 x 10-4. Figure 7 shows the 
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Figure 4: Numerical and exact solution for part a. 


decreasing behavior of the objective function and the rate of convergence in 
the second phase of LP solving. 
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Figure 5: Resulting control function, part b. 


6 Conclusions 


0.8 


0.9 


The method of embedding admissible shapes into a subset of measures is 
extended to an LP formulation for a class of optimal shape design prob- 
lems resulting in a numerical solution scheme. Recent works on this subject 
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Figure 6: Numerical and exact solution for, part b. 
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Figure 7: Objective function in iterations of the simplex method. 


deal with optimization on a specified known part of the shape, while in the 
present paper, the optimization takes place on entirely unknown domain. The 
method is self starting by means that it does not need any initial guess of the 
solution to be started. It is also independent from type of partial differential 
equations and has been presented for a relatively general variational form of 
equations. Numerical examples were used to interpret the implementation of 
the method and to check its validity. 
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